#################################################################################
# Replication file for:                                                         #
# "Balancing Precision and Retention in Experimental Design"                    #
#                                                                               #
# Gustavo Diaz                                                                  #
# Northwestern University                                                       #
# gustavo.diaz@northwestern.edu                                                 #
#                                                                               #
# Erin L. Rossiter                                                              #
# University of Notre Dame                                                      #
# erossite@nd.edu                                                               #
#                                                                               #
# This file produces all tables and figures reported in Appendix Sections       #
# E.3, E.4, E.5, and E.6                                                        #
#################################################################################

# Data ---- 
df_dh <- readRDS("./data/processed_data/DietrichHayesReplication-clean.rds")
df_bg <- readRDS("./data/processed_data/BayramGrahamReplication-clean.rds")
df_th <- readRDS("./data/processed_data/TappinHewittReplication-clean.rds")

# Explicit sample loss (Appendix Table E2) -----

## Design-level sample sizes (column 3)
## Number of units randomized to a design, whether they finished or not
table(df_dh$design, useNA = "ifany")
table(df_bg$design, useNA = "ifany")
table(df_th$design, useNA = "ifany")

## Lost unit at any point (column 4)

### Dietrich & Hayes
table(df_dh$Finished, df_dh$design, useNA = "ifany")
round(prop.table(table(df_dh$Finished, df_dh$design, useNA = "ifany"), margin = 2), 3)

### Bayram and Graham
table(df_bg$Finished, df_bg$design, useNA = "ifany")
round(prop.table(table(df_bg$Finished, df_bg$design, useNA = "ifany"), margin = 2), 3)

### Tappin & Hewitt
### This study is muli-wave, so different.
### TRUE means they attrited somewhere *during* wave 1, 2, or 3
###      or they did not return for wave 2 or 3
df_th$any_loss <- ((df_th$finished_w1 != 1) | 
                   (df_th$finished_w2 != 1) | 
                   (df_th$finished_w3 != 1) | 
                   is.na(df_th$start_w2) | 
                   is.na(df_th$start_w3))
th_loss <- table(df_th$any_loss,
                 df_th$design,
                 useNA = "ifany")
th_loss
round(prop.table(th_loss, margin = 2), 3)


## Specifically post-treatment loss (column 5)

### Dietrich & Hayes
### Because the treatment embedded data was created when the treatment
### screen was shown, if this variable (arbitrarily chosen 1 of 3 treatment indicators)
### is NA, then they attrited *prior* to treatment, if not,
### they attrited post-treatment
tab_dh <- table(df_dh$design, !is.na(df_dh$civilrights), useNA = "ifany")
tab_dh
round(prop.table(tab_dh, margin = 1), 3)



### Bayram and Graham
### Same as DH
tab_bg <- table(df_bg$design, !is.na(df_bg$treatment), useNA = "ifany")
round(prop.table(tab_bg, margin = 1), 3)
tab_bg


### Tappin & Hewitt
### Study is muli-wave, so different.
### Post-treatment attrition is:
### (1) anytime during wave 2 for designs 1-2
### because treatment was immediately assigned when survey started
### (2) anytime during wave 2 or not starting wave 2 for design 3
### because treatment for design 3 was assigned between w1 and w2
### (3) anytime during wave 3
### (4) or not returning to wave 3

### However, some participants were excluded from blocking due to 
### missingness in pre-treatment data.
### These participants definitely count as explicit sample loss,
### but not post-treatment loss. They were not invited back
### to the experimental portion or randomized to treatment.
df_th$pt_attrit <- case_when(df_th$excluded_from_blocking == 1 ~ 0,
                             df_th$finished_w2 == FALSE ~ 1,
                             (df_th$design == 3) & is.na(df_th$start_w2) ~ 1,
                             df_th$finished_w3 == FALSE ~ 1,
                             is.na(df_th$start_w3) ~ 1,
                             .default = 0)    
tab_th <- table(df_th$design, df_th$pt_attrit, useNA = "ifany")
tab_th <- tab_th[,c(2,1)]
tab_th
round(prop.table(tab_th, margin = 1), 3)


# pvalues (column 6)
fisher.test(tab_dh[c(1:2),]) #design 1 vs. 2
fisher.test(tab_dh[c(1,3),]) #design 1 vs. 3

fisher.test(tab_bg[c(1:2),]) #design 1 vs. 2
fisher.test(tab_bg[c(1,3),]) #design 1 vs. 3

fisher.test(tab_th[c(1:2),]) #design 1 vs. 2
fisher.test(tab_th[c(1,3),]) #design 1 vs. 3

tableE2_df <- data.frame(Study = c("Dietrich \\& Hayes", "", "",
                                 "Bayram \\& Graham", "", "",
                                 "Tappin \\& Hewitt", "", ""),
                       Design = c("1. Standard", "2. Pre-post", "3. Pre-post \\& blocking",
                                  "1. Standard", "2. Pre-post", "3. Pre-post \\& blocking",
                                  "1. Standard", "2. Pre-post", "3. Pre-post \\& blocking"),
                       n = c(426, 427, 429,
                             1008, 1010, 1006,
                             751, 780, 767),
                       losses = c("0.005 (2)", "0.007 (3)", "0.007 (3)",
                                  "0.005 (5)", "0.005 (5)", "0.007 (7)",
                                  "0.130 (98)", "0.142 (111)", "0.207 (159)"),
                       pt_loss = c("0.000 (0)", "0.005 (2)", "0.002 (1)",
                                   "0.001 (1)", "0.002 (2)", "0.005 (5)",
                                   "0.125 (94)", "0.141 (110)", "0.132 (101)"),
                       p = c("", "0.499", "1",
                             "", "1", "0.124",
                             "", "0.368", "0.759"))
tableE2 <- kbl(tableE2_df,
    format = "latex",
    align = "l",
    escape = F,
    col.names = c("Study", "Design", "Starting sample size", 'All loss', 'Post-treatment loss', "p-value"),
    booktabs = TRUE,
    linesep = c("","","\\addlinespace"), #for every 3rd row space
    caption = "Explicit Sample Loss by Study and Alternative Design  \\label{tab:exploss}") %>%
  kable_styling(full_width = FALSE, font_size = 10) %>%
  footnote(general = "Descriptive statistics for explicit sample loss across replication studies.
           The starting sample size column shows the number of observations randomized to
           each alternative design. The next two columns show the proportion and count
           in parentheses of how many observations were lost in total and specifically
           post-treatment. The final column shows p-values from Fisher's exact tests of
           the null hypothesis of no difference in proportion comparing post-treatment
           sample loss between each study's standard design and the two alternative designs.",
           footnote_as_chunk = TRUE, threeparttable = TRUE)

writeLines(tableE2, "./tables/tableE2.txt")


# Prep for assessing differential sample inclusion and post-treatment attrition -----

## clean demographics into dummy variables
df_th_mods <- df_th %>%
  mutate(age18_24 = age == 2,
         age25_34 = age == 3,
         age35_44 = age == 4,
         age45_54 = age == 5,
         age55_64 = age == 6,
         age65_74 = age == 7,
         age75_84 = age == 8,
         age85plus = age == 9) %>%
  mutate(gender_m = gender == 1,
         gender_w = gender == 2,
         gender_nb = gender == 3) %>%
  mutate(race_white = grepl(1, race),
         race_black_africanam = grepl(2, race),
         race_amindian_aknative = grepl(3, race),
         race_asian = grepl(4, race),
         race_nativehawaiian_pi = grepl(5, race),
         race_hispanic_latino = grepl(6, race),
         race_arab_middleeastern = grepl(7, race),
         race_other = grepl(8, race)) %>%
  mutate(inc_less20 = inc == 1,
         inc_20_34 = inc == 2,
         inc_35_49 = inc == 3,
         inc_50_74 = inc == 4,
         inc_75_99 = inc == 5,
         inc_100plus = inc == 6) %>%
  mutate(pid_strongdem = ifelse(pid == 1 & pid_strength == 1, 1, 0),
         pid_weakdem = ifelse(pid == 1 & pid_strength == 2, 1, 0),
         pid_leandem = ifelse(pid_lean == 2 & is.na(pid_strength), 1, 0),
         #pid_ind = ifelse(pid_lean == 3, 1, 0),
         pid_leanrep = ifelse(pid_lean == 1 & is.na(pid_strength), 1, 0),
         pid_weakrep = ifelse(pid == 2 & pid_strength == 2, 1, 0),
         pid_strongrep = ifelse(pid == 2 & pid_strength == 1, 1, 0)) %>%
  # political variables asked only in designs 2 and 3
  mutate(ideo_verylib = ifelse(ideo == 1, 1, 0),
         ideo_lib = ifelse(ideo == 2, 1, 0),
         ideo_somewhatlib = ifelse(ideo == 3, 1, 0),
         ideo_moderate = ifelse(ideo == 4, 1, 0),
         ideo_somwhatcons = ifelse(ideo == 5, 1, 0),
         ideo_cons = ifelse(ideo == 6, 1, 0),
         ideo_verycons = ifelse(ideo == 7, 1, 0)) %>%
  rename(ft_voters_reps = ft_voters_pre_1,
         ft_voters_dems = ft_voters_pre_2,
         ft_trump = ft_elites_pre_1,
         ft_obama = ft_elites_pre_2,
         ft_biden = ft_elites_pre_3,
         ft_romney = ft_elites_pre_4)


## pretreatment demographics
## asked in all designs
pretreat_demog <- df_th_mods %>% 
  select(age18_24:pid_strongrep) %>%
  ## remove variables with little to no variation
  select(-age85plus, -gender_nb, -race_other,
         -race_amindian_aknative, -race_nativehawaiian_pi,
         -race_arab_middleeastern) %>%
  colnames()

## pretreatment outcomes & extra political variables
## asked in designs 2 and 3
pretreat_pol <- df_th_mods %>%
  select(ideo_verylib:ideo_verycons,
         ft_voters_reps:ft_romney, 
         salestax_pre,
         pension_pre,
         fedaudit_pre,
         foreignaid_pre,
         healthcare_pre) %>%
  colnames()


# Differential pre-treatment sample inclusion across designs T&H (Appendix Tables E3 and E4) -----

formulas_demog <- lapply(pretreat_demog, function(var) {
  as.formula(paste0(var, " ~ as.factor(design)"))
})

formulas_pol <- lapply(pretreat_pol, function(var) {
  as.formula(paste0(var, " ~ as.factor(design)"))
})

mods_demog <- map(formulas_demog, ~ lm(.x, data = df_th_mods[!df_th_mods$any_loss,]))
mods_pol <- map(formulas_pol, ~ lm(.x, data = df_th_mods[!df_th_mods$any_loss & df_th_mods$design %in% c(2,3),]))

# Because there are too many models, splitting the tables
# into two each, then combining manually as separate rows in the SI
tableE3_top <- stargazer::stargazer(mods_demog[1:13],
                     title = "Assessing Differential Inclusion in the Sample Across Designs - Pre-treatment Demographics Asked in All Designs",
                     label = "tab:diffincldemog",
                     dep.var.caption = "",
                     dep.var.labels.include = F,
                     star.cutoffs = .05,
                     column.sep.width = "0pt",
                     star.char = "*",
                     no.space = T,
                     font.size = "tiny",
                     omit.stat = c("rsq", "adj.rsq", "ser", "f"),
                     covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
                     column.labels =  c("Age (18-24)", "Age (25-34)", "Age (35-44)", "Age (45-54)", "Age (55-64)",
                                       "Age (65-74)", "Age (75-84)", "Gender (Man)", "Gender (Woman)",
                                       "Race (White)", "Race (Black or African American)",
                                       "Race (Asian)", "Race (Hispanic or Latino)"))

writeLines(tableE3_top, "./tables/tableE3_tophalf.txt")


# (second half of table above - use only the tabular part)
tableE3_bottom <- stargazer::stargazer(mods_demog[14:25],
                     title = "",
                     label = "",
                     dep.var.caption = "",
                     dep.var.labels.include = F,
                     star.cutoffs = .05,
                     column.sep.width = "0pt",
                     star.char = "*",
                     no.space = T,
                     font.size = "tiny",
                     omit.stat = c("rsq", "adj.rsq", "ser", "f"),
                     covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
                     column.labels =  c("Inc (<20K)", "Inc (30-34K)", "Inc (35-49K)", "Inc (50-74K)",
                                        "Inc (75-99K)", "Inc (100K+)", "Party (Strong Dem)",
                                        "Party (Weak Dem)", "Party (Lean Dem)",
                                        "Party (Lean Rep)", "Party (Weak Rep)",
                                        "Party (Strong Rep)"))

writeLines(tableE3_bottom, "./tables/tableE3_bottomhalf.txt")

# manually add in row for p-value comparing designs 2 and 3
ftest_pvals <- lapply(mods_demog, car::linearHypothesis, "as.factor(design)2 = as.factor(design)3")
names(ftest_pvals) <- pretreat_demog
ftest_pvals
unlist(lapply(ftest_pvals, function(x) round(x$`Pr(>F)`, 3)[2]))

# Because there are too many models, splitting the tables
# into two each, then combining manually as separate rows in the SI
tableE4_top <- stargazer::stargazer(mods_pol[1:9],
                     title = "Assessing Differential Inclusion in the Sample Across Designs - Additional Pre-treatment Political Questions Asked in Alternative Designs Only",
                     label = "tab:diffinclpol",
                     dep.var.caption = "",
                     dep.var.labels.include = F,
                     star.cutoffs = .05,
                     column.sep.width = "0pt",
                     star.char = "*",
                     no.space = T,
                     font.size = "tiny",
                     omit.stat = c("rsq", "adj.rsq", "ser", "f"),
                     covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
                     column.labels =  c("Ideology (Very Lib.)", "Ideology (Lib.)", "Ideology (Somewhat Lib.)",
                                        "Ideology (Moderate)", "Ideology (Somewhat Cons.)",
                                        "Ideology (Cons.)", "Ideology (Very Cons.)",
                                        "Feeling Thermom. (Voters - Reps)", "Feeling Thermom. (Voters - Dems)"))

writeLines(tableE4_top, "./tables/tableE4_tophalf.txt")

# (second half of table above - use only the tabular part)
tableE4_bottom <- stargazer::stargazer(mods_pol[10:18],
                     title = "",
                     label = "",
                     dep.var.caption = "",
                     dep.var.labels.include = F,
                     star.cutoffs = .05,
                     column.sep.width = "0pt",
                     star.char = "*",
                     no.space = T,
                     font.size = "tiny",
                     omit.stat = c("rsq", "adj.rsq", "ser", "f"),
                     covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
                     column.labels =  c("Feeling Thermom. (Elites - Trump)",
                                        "Feeling Thermom. (Elites - Obama)",
                                        "Feeling Thermom. (Elites - Biden)",
                                        "Feeling Thermom. (Elites - Romney)",
                                        "Salestax Pref.",
                                        "Pension Pref.",
                                        "Fed. Audit Pref.",
                                        "Foreign Aid Pref.",
                                        "Healthcare Pref."))

writeLines(tableE4_bottom, "./tables/tableE4_bottomhalf.txt")


 
# Differential post-treatment attrition across designs T&H (Appendix Figures E1 and E2) -----

## fit models for each pre-treatment covariate and design comparison
## although not identical indicators, removing indicator for woman because few
## gender non-conforming participants, so it is repetitive of "gender_m"
formulas_demog <- lapply(pretreat_demog[-which(pretreat_demog == "gender_w")], function(var) {
  as.formula(paste("pt_attrit ~ as.factor(design) *", var))
})

formulas_pol <- lapply(pretreat_pol, function(var) {
  as.formula(paste("pt_attrit ~ as.factor(design) *", var))
})

mods_demog <- map(formulas_demog, ~ lm_robust(.x, data = df_th_mods))
mods_pol <- map(formulas_pol, ~ lm_robust(.x, data = df_th_mods[df_th_mods$design %in% c(2,3), ]))

## clean and plot for demographic variables 

int_pvalues <- plyr::ldply(.data = mods_demog, .fun = function(x){
  data.frame(p = x$p.value[5:6],
             design = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design"),
             variable = names(x$p)[4], # Variable of interest
             color = if_else(x$coefficients[5:6] > 0, "Larger", "Smaller"),
             row.names = NULL)
})

## No significant interactions -- gender appears significant in
## plot but it is not, reported in Section E.5.1
int_pvalues$p[int_pvalues$variable == "gender_mTRUE"]

# clean for plot
int_pvalues <- int_pvalues %>%
  mutate(design = factor(design, levels = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design")),
         variable_plot = case_when(
           variable == "age18_24TRUE" ~ "Age (18-24)",
           variable == "age25_34TRUE" ~ "Age (25-34)",
           variable == "age35_44TRUE" ~ "Age (35-44)",
           variable == "age45_54TRUE" ~ "Age (45-54)",
           variable == "age55_64TRUE" ~ "Age (55-64)",
           variable == "age65_74TRUE" ~ "Age (65-74)",
           variable == "age75_84TRUE" ~ "Age (74-84)",
           variable == "gender_mTRUE" ~ "Gender (Man)",
           variable == "race_whiteTRUE" ~ "Race/Ethnicity (White)",
           variable == "race_black_africanamTRUE" ~ "Race/Ethnicity (Black or African American)",
           variable == "race_asianTRUE" ~ "Race/Ethnicity (Asian)",
           variable == "race_hispanic_latinoTRUE" ~ "Race/Ethnicity (Hispanic or Latino)",
           variable == "inc_less20TRUE" ~ "Income (<20K)",
           variable == "inc_20_34TRUE" ~ "Income (20-34K)",
           variable == "inc_35_49TRUE" ~ "Income (35-50K)",
           variable == "inc_50_74TRUE" ~ "Income (50-75K)",
           variable == "inc_75_99TRUE" ~ "Income (75-99K)",
           variable == "inc_100plusTRUE" ~ "Income (100K+)",
           variable == "pid_strongdem" ~ "Partisanship (Strong Dem.)",
           variable == "pid_weakdem" ~ "Partisanship (Weak Dem.)",
           variable == "pid_leandem" ~ "Partisanship (Lean Dem.)",
           variable == "pid_leanrep" ~ "Partisanship (Lean Rep.)",
           variable == "pid_weakrep" ~ "Partisanship (Weak Rep.)",
           variable == "pid_strongrep" ~ "Partisanship (Strong Rep.)",
           .default = NA
         ))

## sort based on average p-value per variable for plotting only
avg_pvalues <- as_tibble(int_pvalues) %>%
  group_by(variable_plot) %>%
  summarize(avg_p = mean(p)) %>%
  arrange(avg_p)

int_pvalues <- int_pvalues %>%
  mutate(variable_plot = factor(variable_plot, levels = avg_pvalues$variable_plot))

## plot
plot_attrit <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
  facet_wrap(~design) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
  labs(x = "P-value for the difference in a characteristic's effect on attrition\nbetween the alternative and standard design",
       y = "",
       color = "Effect of characteristic on attrition is\nlarger/smaller in the alternative design") +
  scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(hjust = 1, size = 12),
    panel.border = element_rect(color = "black", fill = NA),
    text = element_text(size = 12),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    strip.text = element_text(size = 12),
    legend.position = "bottom"
  ) +
  scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
                     labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

## save
ggsave("figures/figE1.pdf", plot_attrit, width = 10, height = 6)


## clean and plot for political variables 

int_pvalues <- plyr::ldply(.data = mods_pol, .fun = function(x){
  data.frame(p = x$p.value[4],
             design = c("Pre-post & blocking vs.\nstandard design"),
             variable = names(x$p)[3], # Variable of interest
             color = if_else(x$coefficients[4] > 0, "Larger", "Smaller"),
             row.names = NULL)
})


# clean for plot
int_pvalues <- int_pvalues %>%
  mutate(variable_plot = case_when(
           variable == "ideo_verylib" ~ "Ideology (Very Liberal)",
           variable == "ideo_lib" ~ "Ideology (Liberal)",
           variable == "ideo_somewhatlib" ~ "Ideology (Somewhat Liberal)",
           variable == "ideo_moderate" ~ "Ideology (Moderate)",
           variable == "ideo_somwhatcons" ~ "Ideology (Somewhat Conservative)",
           variable == "ideo_cons" ~ "Ideology (Conservative)",
           variable == "ideo_verycons" ~ "Ideology (Very Conservative)",
           variable == "ft_voters_reps" ~ "Feeling Thermometer (Voters - Reps)",
           variable == "ft_voters_dems" ~ "Feeling Thermometer (Voters - Dems)",
           variable == "ft_trump" ~ "Feeling Thermometer (Elites - Trump)",
           variable == "ft_obama" ~ "Feeling Thermometer (Elites - Obama)",
           variable == "ft_biden" ~ "Feeling Thermometer (Elites - Biden)",
           variable == "ft_romney" ~ "Feeling Thermometer (Elites - Romney)",
           variable == "salestax_pre" ~ "Salestax Preference",
           variable == "pension_pre" ~ "Pension Preference",
           variable == "fedaudit_pre" ~ "Fed. Audit Preference",
           variable == "foreignaid_pre" ~ "Foreign Aid Preference",
           variable == "healthcare_pre" ~ "Healthcare Preference",
           .default = NA
         ))

## sort based on p-value size
sort_pvalues <- as_tibble(int_pvalues) %>%
  group_by(variable_plot) %>%
  arrange(p)

int_pvalues <- int_pvalues %>%
  mutate(variable_plot = factor(variable_plot, levels = sort_pvalues$variable_plot))


## plot
plot_attrit_pol <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
  labs(x = "P-value for the difference in a characteristic's effect on attrition\nbetween the two alternative designs",
       y = "",
       color = "Effect of characteristic on attrition is larger/smaller\nin pre-post with blocking vs. pre-post alone") +
  scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(hjust = 1, size = 12),
    panel.border = element_rect(color = "black", fill = NA),
    text = element_text(size = 12),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    strip.text = element_text(size = 12),
    legend.position = "bottom"
  ) +
  scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
                     labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

## save
ggsave("figures/figE2.pdf", plot_attrit_pol, width = 10, height = 6)



# Differential post-treatment attrition across treatment arms by design T&H (Appendix Figures E3 and E4) -----

## fit models for each pre-treatment covariate and design comparison
## although not identical indicators, removing indicator for woman because few
## gender non-conforming participants, so it is repetitive of "gender_m"

## also no variation in Age 74-84
formulas_demog <- lapply(pretreat_demog[-which(pretreat_demog %in% c("gender_w", "age75_84"))], function(var) {
  as.formula(paste("pt_attrit ~ as.factor(design) * Z *", var))
})

formulas_pol <- lapply(pretreat_pol, function(var) {
  as.formula(paste("pt_attrit ~ as.factor(design) * Z *", var))
})

mods_demog <- map(formulas_demog, ~ lm_robust(.x, data = df_th_mods))
mods_pol <- map(formulas_pol, ~ lm_robust(.x, data = df_th_mods[df_th_mods$design %in% c(2,3), ]))

## clean and plot for demographic variables 

int_pvalues <- plyr::ldply(.data = mods_demog, .fun = function(x){
  data.frame(p = x$p.value[11:12],
             design = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design"),
             variable = names(x$p)[5], # Variable of interest
             color = if_else(x$coefficients[11:12] > 0, "Larger", "Smaller"),
             row.names = NULL)
})


# clean for plot
int_pvalues <- int_pvalues %>%
  mutate(design = factor(design, levels = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design")),
         variable_plot = case_when(
           variable == "age18_24TRUE" ~ "Age (18-24)",
           variable == "age25_34TRUE" ~ "Age (25-34)",
           variable == "age35_44TRUE" ~ "Age (35-44)",
           variable == "age45_54TRUE" ~ "Age (45-54)",
           variable == "age55_64TRUE" ~ "Age (55-64)",
           variable == "age65_74TRUE" ~ "Age (65-74)",
           variable == "age75_84TRUE" ~ "Age (74-84)",
           variable == "gender_mTRUE" ~ "Gender (Man)",
           variable == "race_whiteTRUE" ~ "Race/Ethnicity (White)",
           variable == "race_black_africanamTRUE" ~ "Race/Ethnicity (Black or African American)",
           variable == "race_asianTRUE" ~ "Race/Ethnicity (Asian)",
           variable == "race_hispanic_latinoTRUE" ~ "Race/Ethnicity (Hispanic or Latino)",
           variable == "inc_less20TRUE" ~ "Income (<20K)",
           variable == "inc_20_34TRUE" ~ "Income (20-34K)",
           variable == "inc_35_49TRUE" ~ "Income (35-50K)",
           variable == "inc_50_74TRUE" ~ "Income (50-75K)",
           variable == "inc_75_99TRUE" ~ "Income (75-99K)",
           variable == "inc_100plusTRUE" ~ "Income (100K+)",
           variable == "pid_strongdem" ~ "Partisanship (Strong Dem.)",
           variable == "pid_weakdem" ~ "Partisanship (Weak Dem.)",
           variable == "pid_leandem" ~ "Partisanship (Lean Dem.)",
           variable == "pid_leanrep" ~ "Partisanship (Lean Rep.)",
           variable == "pid_weakrep" ~ "Partisanship (Weak Rep.)",
           variable == "pid_strongrep" ~ "Partisanship (Strong Rep.)",
           .default = NA
         ))

## sort based on average p-value per variable for plotting only
avg_pvalues <- as_tibble(int_pvalues) %>%
  group_by(variable_plot) %>%
  summarize(avg_p = mean(p)) %>%
  arrange(avg_p)

int_pvalues <- int_pvalues %>%
  mutate(variable_plot = factor(variable_plot, levels = avg_pvalues$variable_plot))

## plot
plot_attritZ <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
  facet_wrap(~design) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
  labs(x = "P-value for the difference in a characteristic's effect\non differential attrition across treatment arms\nbetween the alternative and standard design",
       y = "",
       color = "Effect of characteristic on attrition is\nlarger/smaller in the alternative design") +
  scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(hjust = 1, size = 12),
    panel.border = element_rect(color = "black", fill = NA),
    text = element_text(size = 12),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    strip.text = element_text(size = 12),
    legend.position = "bottom"
  ) +
  scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
                     labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

## save
ggsave("figures/figE3.pdf", plot_attritZ, width = 10, height = 6)


## clean and plot for political variables

int_pvalues <- plyr::ldply(.data = mods_pol, .fun = function(x){
  data.frame(p = x$p.value[8],
             design = c("Pre-post & blocking vs.\nstandard design"),
             variable = names(x$p)[4], # Variable of interest
             color = if_else(x$coefficients[8] > 0, "Larger", "Smaller"),
             row.names = NULL)
})


# clean for plot
int_pvalues <- int_pvalues %>%
  mutate(variable_plot = case_when(
    variable == "ideo_verylib" ~ "Ideology (Very Liberal)",
    variable == "ideo_lib" ~ "Ideology (Liberal)",
    variable == "ideo_somewhatlib" ~ "Ideology (Somewhat Liberal)",
    variable == "ideo_moderate" ~ "Ideology (Moderate)",
    variable == "ideo_somwhatcons" ~ "Ideology (Somewhat Conservative)",
    variable == "ideo_cons" ~ "Ideology (Conservative)",
    variable == "ideo_verycons" ~ "Ideology (Very Conservative)",
    variable == "ft_voters_reps" ~ "Feeling Thermometer (Voters - Reps)",
    variable == "ft_voters_dems" ~ "Feeling Thermometer (Voters - Dems)",
    variable == "ft_trump" ~ "Feeling Thermometer (Elites - Trump)",
    variable == "ft_obama" ~ "Feeling Thermometer (Elites - Obama)",
    variable == "ft_biden" ~ "Feeling Thermometer (Elites - Biden)",
    variable == "ft_romney" ~ "Feeling Thermometer (Elites - Romney)",
    variable == "salestax_pre" ~ "Salestax Preference",
    variable == "pension_pre" ~ "Pension Preference",
    variable == "fedaudit_pre" ~ "Fed. Audit Preference",
    variable == "foreignaid_pre" ~ "Foreign Aid Preference",
    variable == "healthcare_pre" ~ "Healthcare Preference",
    .default = NA
  ))

## sort based on p-value size
sort_pvalues <- as_tibble(int_pvalues) %>%
  group_by(variable_plot) %>%
  arrange(p)

int_pvalues <- int_pvalues %>%
  mutate(variable_plot = factor(variable_plot, levels = sort_pvalues$variable_plot))


## plot
plot_attritZ_pol <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
  labs(x = "P-value for the difference in a characteristic's effect\non differential attrition across treatment arms\nbetween the two alternative designs",
       y = "",
       color = "Effect of characteristic on attrition is larger/smaller\nin pre-post with blocking vs. pre-post alone") +
  scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(hjust = 1, size = 12),
    panel.border = element_rect(color = "black", fill = NA),
    text = element_text(size = 12),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    strip.text = element_text(size = 12),
    legend.position = "bottom"
  ) +
  scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
                     labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

## save
ggsave("figures/figE4.pdf", plot_attritZ_pol, width = 10, height = 6)


# Treatment effect estimates with full samples (Appendix Tables E5, E6, E7) -----

## Design 1
dh_design1 <- lm_robust(speech_approve_scaled ~ Z, data = df_dh, subset = design == 1)

bg_design1 <- lm_robust(delpref_post ~ Z, data = df_bg, subset = design == 1)

th_design1 <- lm_robust(pension_followup_scaled ~ Z, data = df_th, subset = design == 1)


## Design 2
dh_design2 <- lm_robust(speech_approve_scaled ~ Z + dv_quasipre_scaled, data = df_dh, subset = design == 2)

bg_design2 <- lm_robust(delpref_post ~ Z + delpref_pre, data = df_bg, subset = design == 2)

th_design2 <- lm_robust(pension_followup_scaled ~ Z + pension_pre_scaled, data = df_th, subset = design == 2)


## Design 3
dh_design3 <- lm_robust(speech_approve_scaled ~ Z, fixed_effects = ~ block, data = df_dh, subset = design == 3)

bg_design3 <- lm_robust(delpref_post ~ Z, fixed_effects = ~ block, data = df_bg, subset = design == 3)

th_design3 <- lm_robust(pension_followup_scaled ~ Z + pension_pre_scaled, fixed_effects = ~ block_id, data = df_th, subset = design == 3)

## Pooled treatment effects
dh_pooled <- lm_robust(speech_approve_scaled ~ Z, data = df_dh)

bg_pooled <- lm_robust(delpref_post ~ Z, data = df_bg)

th_pooled <- lm_robust(pension_followup_scaled ~ Z, data = df_th)


# Create tables

dh_mod_list <- list("Pooled" = dh_pooled, "Standard design" = dh_design1, "Pre-post" = dh_design2, "Block randomized" = dh_design3)

bg_mod_list <- list("Pooled" = bg_pooled, "Standard design" = bg_design1, "Pre-post" = bg_design2, "Block randomized" = bg_design3)

th_mod_list <- list("Pooled" = th_pooled, "Standard design" = th_design1, "Pre-post" = th_design2, "Block randomized" = th_design3)

tableE5 <- modelsummary::modelsummary(dh_mod_list,
                           coef_map = c("Z" = "Treatment Effect",
                                        "dv_quasipre_scaled" = "Pre-treatment outcome measure",
                                        "(Intercept)" = "Intercept"),
                           output = "latex",
                           stars = c('*' = .05),
                           gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0)),
                           add_rows = data.frame("Block fixed effects", "", "", "", "Yes"),
                           title = "Treatment Effect Estimates By Design, Dietrich and Hayes Replication \\label{tab:dhests}",
                           escape = F)

tableE6 <- modelsummary::modelsummary(bg_mod_list,
                           coef_map = c("Z" = "Treatment Effect",
                                        "delpref_pre" = "Pre-treatment outcome measure",
                                        "(Intercept)" = "Intercept"),
                           output = "latex",
                           stars = c('*' = .05),
                           gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0)),
                           add_rows = data.frame("Block fixed effects", "", "", "", "Yes"),
                           title = "Treatment Effect Estimates By Design, Bayram and Graham Replication \\label{tab:bgests}",
                           escape = F)

tableE7 <- modelsummary::modelsummary(th_mod_list,
                           coef_map = c("Z" = "Treatment Effect",
                                        "pension_pre_scaled" = "Pre-treatment outcome measure",
                                        "(Intercept)" = "Intercept"),
                           output = "latex",
                           stars = c('*' = .05),
                           gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0)),
                           add_rows = data.frame("Block fixed effects", "", "", "", "Yes"),
                           title = "Treatment Effect Estimates By Design, Tappin and Hewitt Replication \\label{tab:thests}",
                           escape = F)

writeLines(as.character(tableE5), "./tables/tableE5.txt")
writeLines(as.character(tableE6), "./tables/tableE6.txt")
writeLines(as.character(tableE7), "./tables/tableE7.txt")


